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Abstract 

A substantial increase of the mean logarithmic mass (In A) of galactic cosmic rays vs 
energy has been observed . We study three effects that could explain this trend i) dif- 
ferent source spectra for protons and heavy nuclei ii) a selective nuclear destruction 
in flight of heavies iii) a gradient of the source number and chemical composition 
in the galactic disk. We take advantage of the diffusive cosmic ray propagation 
model developed at LAPTH to study specifically the geometrical aspects of the 
propagation and extend it to high energy. Using a simple modeling of the spectral 
knee around 10 15 eV, a bump in (In A) appears. This feature is smoother when the 
spectral index of protons is steeper than Fe's. We analyze the effects of the rigidity 
dependence of the diffusion coefficient and the scale height of the confinement halo 
and we show that (hi A) is most sensitive to the first parameter. Pure geometrical ef- 
fects are less determining than the diffusion coefficient spectral index. Subsequently, 
we conclude that the physics of cosmic ray confinement is the essential cause of the 
heavy nuclei enrichment until ~ 10 15 eV. 
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Introduction 



Recent measurements of the cosmic ray average logarithmic mass and all- 
particle spectrum around 10 15 eV [1,2] give new clues to understand the origin 
of the cosmic rays and in particular the puzzle of the knee in the energy spec- 
trum. Analysis of such data in a coherent theoretical framework is a rough 
task, and even if much progress has been done both in theoretical and exper- 
imental sides, none of the models proposed so far to solve this problem has 
been unanimously accepted. The highest energy particles are almost certainly 
extragalactic. A similar origin is not excluded near the knee, but it is difficult 
to account for the observed continuity of the spectrum in this region [3]. As 
a consequence, the intermediate region between 10 15 and 10 19 eV should be 
analysed in terms of the same physical mechanisms than lower energy parti- 
cles. 

There are at least three explanations for the knee: (i) a change in propagation 
parameters (diffusive regime), (ii) a change in the source regime, (iii) a change 
related to the properties of high energy interactions in the atmosphere or a 
subtle combination of all three. As Schatz [4] emphasized recently, the fine 
structure of the knee in all-particle spectra provided by extensive air showers 
can help to discriminate between these solutions. Furthermore, recent data 
from collider seem to show no drastic departure from cross sections predictions 
in the range 100 - 1000 TeV [1]. Anyway, this paper will only concentrate on 
astrophysical aspects. 

More information can be obtained by measuring (In A) . All-particle spectrum 
and (In A) are given by linear combination of the individual fluxes with dif- 
ferent weights; 
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and 



(ha. A) 



Therefore they provide different information of these fluxes. Obviously, these 
weighted quantities are not very useful at "low energy" where all nuclei are well 
resolved in satellite or balloon experiments (for a compilation of data, see [5]). 
Experimental difficulties arise around and above the knee: fluxes are very low 
(~ m -2 sr" 1 yr _1 ) and large ground array detectors are needed to collect 
unresolved events with a good statistic. Incidently, the all-particle spectrum 
can be extracted via shower parameters (e.g. core position, direction,...) and 
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an estimation of the average logarithmic mass is also possible through various 
methods (e.g [2,6,7] and in particular [8]). A new experimental technique has 
been recently proposed, which can potentially yield excellent charge resolution 
measurements near and above the knee [9]. But at the present time, even if 
a few solutions exist to infer the composition [10,11], all-particle flux, mean 
logarithmic mass and sometimes proton and helium spectra [12] are grossly 
the only available observables near and above the knee (for a detailed revue, 
see [13]). 

The paper is organized as follows: (i) features of primary species propagated 
in diffusion model are reminded; (ii) simple models are used to explore the be- 
haviour of (In A) considering separately the effect of three parameters (source 
spectra, propagation, geometrical aspects); (iii) the results of these simplified 
models are analysed; (iv) using a more realistic simulation, the three parame- 
ters are discussed in details, and (v) all these parameters are included in the 
propagation model of Maurin et al. [14] and are finally combined to analyze 
the knee problem. 



1 Basic features of propagation models 

For practical reasons, spectra at low energy are almost always displayed in 
units of kinetic energy per nucleon, because this quantity is conserved in nu- 
clear reactions. This convention has become the rule in the analysis of cosmic 
ray propagation, e.g. secondary to primary ratio studies. Above PeV ener- 
gies, fluxes are plotted vs energy per particle because it is the only observable 
provided by ground observatories. Whereas various primary species present 
similar spectra when plotted versus kinetic energy per nucleon or rigidity, 
they behave differently in terms of total energy (compare for example Fig. 5-a 
of [15] that displays fluxes in units of kinetic energy per nucleon and Fig. 2 
of [5] that shows the same fluxes, but in units of total energy); so a combina- 
tion of all the corresponding spectra into a single quantity (for example the 
average logarithmic mass) will provide a different dependence in terms of total 
energy or in terms of rigidity. 

1.1 Propagation models above the knee and extrapolation to higher energies 

Although cosmic ray properties are well understood up to a few hundreds 
of GeV/nuc (see the recent review of [16]), it is difficult to harmonize all 
the observables {i.e. protons, nuclei, e~, e + , p, 7 rays; see for instance [17]). 
In particular, the hypothesis according to which the fluxes measured locally 
are representative of the fluxes present everywhere in the Galaxy is still de- 
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bated. At intermediate energies, extension of the usual propagation models 
has mainly to face the problem that the observed anisotropy is very low. This 
seems to favor a Kolmogorov spectrum (interstellar turbulence) associated to 
the rigidity dependence of the diffusion coefficient (see discussion in [18]). 

Cosmic ray spectra are affected by the propagation process. Maurin et al. [14] 
have recently developed a propagation model (see Sec. 4.2 for details) that 
is in principle valid for a wide energetic range as long as charged nuclei are 
considered. This is a two zones model (thin gazeous disc - half-height h = 100 
pc - and large diffusive halo L ~ 3 — 10 kpc) where cylindrical symmetry is 
assumed (radial extension, R = 20 kpc). The steady-state differential density 
Ni(E, r) of the nucleus j as a function of energy E and position (r, z) in this 
model is given by 



L dlff N 3 (r, z) + 2hS(z) ( q 3 Q J (E)q(r) + £ T kj N k (r, 0) - FA^r, 0) 

fc=i 

2hS(z) A lv{E)N J {0) - d\E)-^NM 



d ( d 2 I d d \ 

Ldiff ^ c Qz ^ \dz 2 r dr dr J 

The operator Ldiff represents convection (V c ) plus spatial diffusion (K(E) = 
Kq (3xR 5 ) acting in the whole box. All other terms describe processes localized 
in the thin gazeous disc only: the bracket corresponds to primary source term, 
secondary spallative production from heavier nuclei, and destruction cross 
section (radioactive-induced processes have been omitted here) . Curly bracket 
provides all terms leading to energetic redistribution, i.e losses (Coulombian, 
ionisation and adiabatic expension losses) and gains (reacceleration) described 
by two effective parameters bj(E) and d?{E) (for all details, see [14]). 

In the leaky box model, widely used, all the quantities are spatially averaged. 
The diffusion-convection term is then replaced by an effective escape term 
which has the meaning of a residence time r esc (~ 20 Myr at a few GeV/nuc) 
in the confinement volume (—Ldiff ^ T esl)- ^ energy gains and losses are 
discarded, the leaky box equation reads 

+ UqAe) + E r k >N k - pnA = o 

Tesc y k=1 J 

As an immediate consequence, leaky box models do not allow to take in con- 
sideration, for example, a radial dependence of Galactic sources. More gen- 
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erally, any subtle effect related with spatial dependence of any terms of the 
diffusion-convection equation is automatically swept away. 




10 1 1 10 10 2 io 3 10 4 io : 



E k/nuc [GeV/ nuc] 

Fig. 1. Unnormalized Fe flux for a typical propagation model that reproduces B/C 
(reference curve); the five parameters (diffusion - normalization Kq and slope 5 -, 
halo size L, convection V c , reacceleration V a ) are taken from [14]. Other curves 
correspond to the same propagation parameters, except one that has been switched 
off. 

Figure 1 displays the iron flux as a function of kinetic energy per nucleon in our 
diffusion-convection model; this behaviour is standard for all charged nuclei 
except protons, for which propagation is less sensitive to nuclear reactions. 
Within realistic propagation models, several parameters affect the primeval 
power law source spectrum: 

• The most spectacular effect is that of inelastic interactions (nuclear de- 
struction). It becomes negligible (< a few %) between 1 TeV/nuc and 
10 TeV/nuc, depending on the species, since the escape time is shorter 
than the nuclear time. 

• The spectrum is also affected by energy losses up to ~ 100 GeV/nuc. 

• More controversial effects may be present: galactic wind up to 1 TeV/nuc 
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and reacceleration up to some tens of GeV/nuc. 

Propagation effects are present up to fO TeV/nuc. This induces modifications 
of spectra up of a few hundreds TeV in total energy, at least for heavier nuclei 
(i.e. Fe) leading in turn to a dependence on energy of the average logarithmic 
mass till these higher energies, i.e. near the knee (~ PeV). Since current high 
energy data are provided in total energy per particle, all results will be given 
in total energy (denoted E) unless stated otherwise. 



2 Separation of key ingredients 



Solar modulation is ignored because it has almost no effect on fluxes above a 
few tens of GeV/nuc. Three items (source spectra, propagation and geomet- 
rical effects) are pieces of the cosmic ray puzzle. In this section, we introduce 
simplified models where only two species are considered (e.g. H and Fe) that 
allow to analyze qualitatively the increase of the mean mass of cosmic rays 
with energy. 



2.1 Source spectrum effect 



In the simplest propagation model that one can imagine, pure diffusion is 
considered and only two ingredients are necessary: 

Qj(E) — $ x I ) is the source spectrum, (1) 

VI GV/ 

R ^ 



K(E) = K x ( — — — ) is the diffusion coefficient. (2) 
. 1 G V / 



In the above expressions, R = p/Z is the rigidity, g° and <x,- are respectively the 
source abundance and the spectral index of nucleus j, K is the normalization 
coefficient of the diffusion coefficient and 5 its slope (both are assumed inde- 
pendent of the species considered). The detected local flux (in arbitrary units) 
is thus expressed merely as the ratio of sources over the diffusion coefficient 
(up to a coefficient v/An ~ c/Att for all nuclei): 

$ J (E)=Q J (E)/K(E). (3) 



This model can be viewed as the high energy limit of all diffusion models. This 
limit is also that of a locally equivalent leaky box model, where once emitted, 
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the nucleus can escape from the confinement box with a probability which is 
related very simply to K(E) - see for example [14] and references therein. 



2.1.1 (In A) for the two nuclear component model 

In a two nuclei model, the evaluation of the Fe/p ratio and the average log- 
arithmic mass number - "average mass" for short - is straightforward. At a 
given energy, the rigidities of two charged nuclei are different. Indeed, the term 
K(E) can't be simply factorized, and some residual multiplicative constant de- 
pending on 5 appears. Using the approximate relation R = p/Z ~ EjZ (E is 
the total energy per particle) which is correct beyond a few tens of A GeV, 
we have 

( ya Fc +5\ / W \ - Q Fc+a P 



where the first term corresponds to the relative source abundance Fe/p taken 
at the energy 100 GeV, which is ~ 1/20 (see Fig. 2 of compilation of [5]). 
Apart from numerical constant, we roughly see that the steeper — a Fo + a P 
effective slope, the greater will be the average mass. To be more precise, 



(\nA)(E) 



\n(A 



Fc J 



100 GeV 



x (E/100 GeV)- a p+ a ^ 



"(4) 



whereas in terms of the all-particle spectrum, the evolution - up to a normal- 
ization - can be written as 



$a ii K E -a Fc 



„0 7« p +(5 

„0 yctF c +5 
<?Fc Z : 



Fc 



X 



J 100 GeV 



E 



100 GeV 



-a p +QFe 



• (5) 



We remark that in our simple two species model, $ al1 oc E aFc 5 /(In A)(E); 
the behaviour of all-particles spectrum will be discussed in Sec. 3. 



2.1.2 Application to the"wind/ISM" supernova model 

Among the various acceleration models, some of them are able to produce 
various spectra for different nuclei (see Sec. 4.1). The maximal source ef- 
fect is obviously obtained when | — a p + ap e \ is maximal. For typical val- 
ues — a p + ape ~ —0.1 such as those advocated by the wind/ISM super- 
novae model of Biermann and collaborators [19], and with the above value 
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~ 20, we obtain 

lOOGcV ' 



(q° P /q° Fc )(z^ +5 /z^ +s ] 

(lnA)(E = 100 GeV/ 10 TeV/ 1 PeV/ oo) ~ 0.19 / 0.30 / 0.50 / 4.03 (6) 



Indeed, the evolution is expected to be smoother when all cosmic ray nuclei 
are considered, because the slopes of nuclei heavier than hydrogen are similar 
(see Sec. 4.1); but this sole effect leads to an evolution of the average mass. 
Incidently, if a p = a Fe , we are left with a constant (In -4) that merely depends 
on the relative abundances of nuclei. This conclusion holds again when several 
nuclei are considered. 



2.2 Propagation effect: same spectral index 



The process of inelastic collisions - which differs from one nucleus to another 
- has a major impact on propagation (see Fig. 1). To obtain order of magni- 
tude estimates, a simple leaky box description including just destruction cross 
sections is sufficient. Thus, the propagated flux for a primary species j is: 

^ (n _ QAE) 1 Qj(E) 

A) l/\ esc {E) + af d /A A o- esc {E) + o-T el 



Here aj 1 ^ = a\ ot — of is the total inelastic (or reaction) cross section of species 
j, and X esc (E) is the usual escape length of the leaky box in g cm -2 and A 
is the mean mass of the atoms in the interstellar medium. This can also be 
rewritten in terms of an effective escape cross section a esc (in cm 2 ) that is 
related to the usual grammage by a simple factor. 

Under the assumption of the same source spectrum (i. e. universal slope <x, = a 
for all species) and since there is no geometry in the leaky box, assuming an 
escape dependence X esc (E) = \ fiR~ 5 ~ \ R~ S , i.e a esc (E) = o~oR s for all 
nuclei, we get: 



(lnA)(E) = \n(A Fc )x 



1 + We) 



Z p \ a (a (E/Z Fe ) s + a™ 1 )^^ 



x 



F J ' (a (E/Z p ) s + a 



(8) 



Contrary to the previous case, the ratio q p /qp e does not depend on energy. 
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2.2.1 Leaky box results 

The typical diffusion coefficient for leaky box models is given for example 
by [20] 

A = 35.1/5 (R/l GV)-°' 61 g cm" 2 

which corresponds in the above expression to <j ~ 47 mb. We have in this 
model a ~ 2.2, a™ el ~ 30 mb (Particle Data GroupQ), a l £ el ~ 710 mb [21]. 
In order to compare this effect with the precedent, we set q^/q^c to ~ 25000 
which gives the same (In A) at 100 GeV as before. Thus, 

(\nA){E = 100 GeV/ 10 TeV/ 1 PeV / oo) ~ 0.19 / 0.83 / 1.05 / 1.08 (9) 



Even near the knee energies, the average mass slightly evolves due to spalla- 
tions. It is then clear from formula (8) that the evolution of chemical abun- 
dances depends on (i) the destruction cross section, (ii) the diffusion coefficient 
slope 5, as well as other propagation parameters (see Sec. 1.1). 



2.3 Geometrical effects 



In homogeneous diffusion models {e.g. leaky box models), this would be the 
final step. In the context of diffusion models, the nuclei propagate in a two 
zones/three dimensional space, and we take into account the gradients of the 
source number and of the metallicity. 



2.3.1 Basic description 

Let first assume that spectral source indexes are similar for all accelerated 
species, and that pure diffusion prevails. The quantity 

MW= TTt¥^ > (10) 

1 + (<?p/<?Fe) 



depends on energy via the relative abundance (q p /qFe) ca = (<7p/<7Fb)- Actually, 
the source term depends on radial coordinate, and the average distance (rx) 
from which a given species X come depends on energy. At a given energy, 
(rx) is generally different for two species, implying an indirect dependence 
of (q P /qFe) eS = (o'p/o'Fc) on energy. Moreover, even if these average distances 

2 http://pdg.lbl.gov/ 
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(r p )(E) and (r-F C )(E) are equal, the relative abundances may depend on this 
distance; this is the pure metallicity effect. At sufficiently high energy, i.e. a 
few tens of TeV/nuc, (rx) is a constant number that only depends on the size 
of the diffusive box: no geometrical effects are expected. 



2.3.2 Metallicity effect 

We focus on the pure metallicity effect: forgetting for a while the radial source 
distribution effect, we suppose that relative distribution of species only de- 
pends of their location. Namely, we take a gradient V[Fe/H] ~ —0.05 dex 
kpc" 1 , corresponding to an increased metallicity towards the center of the 
Galaxy (see Sec. 4.4). The sun is located at r = R & (8.5 kpc). The metallicity 
gradient is given by: 

[Fe/H] = -0.05 (r - R & ) = log (q Fc /q p ) cG - log(Fe/H) 



leading to 



c , 







A rough estimation of the r dependent term above can be obtained as follows: 
at high energy, cosmic ray nuclei cannot come farther than the center, i.e. 
r = kpc. At lower energy (say 100 GeV), they all come from r = R@. It 
corresponds to 

0.32 




This crude evaluation that overestimates some effects gives us the evolution of 
the mean mass for the metallicity gradient effect: using Eq. (10) and assuming 
once again that at 100 GeV, the ratio (q p /qFe) cS ~ 20, we obtain 

(\aA)(E = 100 GeV/ oo) ~ 0.19 / 0.54 (12) 



This pure metallicity effect is small compared to the others. Moreover, as 
cosmic ray sources are distributed in the anticenter (lower metallicity) as well 
as in the center (higher metallicity), we expect a sort of cancelation when both 
contributions are added. The same balance is likely to occur for the second 
geometrical effect. 
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3 First conclusions 



Several qualitative remarks can be made about the relative importance of 
the various effects. The metallicity gradient plays an almost negligible role 
compared to the two others. Pure propagation effect is dominant (see Fig 2, 
left panel; right panel shows the two separate effects for the all particle flux), 
but as well as geometrical effects it ceases to act around 100 TeV - 1 PeV 
(roughly the knee energy). Then, the pure source effect only, i.e. different a 
for protons and other species (a p ^ ctz>i), is able to produce an evolution 
of the average mass. This is an important result that validates the approach 
used in a recent study [2]. 




Fig. 2. Two nuclei model, p and Fe. Left panel: (In A) for pure spectrum effect - 
solid line (see Sec. 2.1.2) and pure propagation effects - dashed line (see Sec. 2.2.1). 
Right panel: same for all-particle flux (arbitrary units). 

Including all the primaries, we should expect concerning the first effect - pure 
source spectrum, see Eq. (6) - a smoother evolution, since other primaries show 
slopes similar to that of iron and mix their effects. Concerning the second 
effect, a smoothing is also expected since all nuclei are equally dispatched 
between proton and iron (from the cross section point of view) . 

To visualise the changes induced by the knee, we display in Fig. 3 (for the 
two previous cases) three forms for the break in spectra, namely a break at 
10 15 eV either in total energy, or rigidity or energy per nucleus. Concerning 
(In A), we see that the evolution from lighter to heavier nuclei at the knee 
is more important if the slope of hydrogen and other nuclei are the same 
below the knee. If not the case, we see in particular that the farthest the 
break occurs, the smoother is the bump in (In A) (compare dot line - -Eknee = 
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10 14 eV - to solid line - E^ ncc = 10 15 eV). For each model, we have after the 
knee a constant composition related to the fact that all nuclei have now the 
same slope. Concerning the all-particle spectrum, we see that, except for the 
situation where the break is energy dependent, there is a smooth evolution 
on about one decade before reaching the definitive slope. A first break occurs 
when protons change slope, and the second when Fe does. The situations where 
.Eknee ^ A eV or -Eknee oc Z eV are quite similar, and the difference will be 
rather difficult to see in data. Note that a change of diffusive regime (with 
Ad = 0.4) is completely equivalent to the case Aa = 0.4 in source spectra if 
-E'knee Z eV. This is because at these energies, subtleties of propagation are 
irrelevant. 




-2 
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Fig. 3. Left panel: (In A), and right panel: all-particle flux times E 2 ' 5 (arbitrary 
units). The various cases correspond to the two effects previously studied (see Fig. 2) 
with an additional modeling for the knee. For each curve (dashed and solid line) 
we have from bottom to top a transition £/k nee = 10^ eV, Ek Qee = Z x 10 15 eV, 
and -Ekncc = A x 10 15 eV. The dot curve in left panel corresponds to the case 
£ knee = Zx 10 14 eV. 



4 Cosmic ray diffusion model 

4-1 Spectral indexes of sources 
4-1.1 Acceleration models 

The maximal energy reached in shocks associated with supernovae is of about 
E ~ Z x 10 14 eV (see for instance [22]). As the acceleration processes are 
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rigidity dependent, there is a cut-off at Z times the maximum energy gained 
and these models predict an increase of the average mass of primary cosmic 
rays with total energy. In fact, the limit E ~ Z x 10 14 eV is close to the knee 
energy so that another acceleration process must be found for higher energies. 
To bypass this limit and explain slopes near the knee, several mechanisms 
have been proposed ; for example, acceleration by terminal shock of galactic 
wind [23,24], by neutron star quakes [25] or by pulsars [26], contribution from 
a single recent local supernova explosion [27,28], and photodisintegration of 
nuclei by a background of optical and soft UV photons in the sources [29]. 
Another explanation also has been proposed, based on the extension of super- 
nova acceleration models, but they are contradictory: it is shown in [30] that 
acceleration by multiple spherical shocks in OB associations (superbubbles) 
is not sufficient to reach E ~ 10 17 — 10 18 eV unless extreme values of the 
turbulence parameters are used. A different conclusion is drawn in [31], where 
the authors adjust turbulence parameters to reproduce the data. Recently, 
going back to the problem of the maximal energy gained in SN shocks, [32,33] 
showed that non-linear amplifications of the magnetic field by the cosmic rays 
themselves could push the usual limit to Z x 10 17 eV, and even a factor ten 
more if stellar wind pre-exists. Note that this distinction between standard 
supernovae and supernova explosions of massive stars in their own wind has 
been advocated by Biermann and collaborators [34,35] as the possible expla- 
nation of the knee. This a is very important argument because it seems that 
such models produce below the knee different spectra for p and other species. 
Note that even in usual acceleration models, collective effects can also produce 
such an effect [36]. 

4-1.2 Data: behaviour at the knee 

The energy of the break and slopes below and above the knee (denoted ji and 
72) are in relative agreement between the various experiments (see Tab. 1). 
Direct measurements of fluxes and extraction from air showers experiments of 
proton and helium fluxes give a constant slope at least till Zx 10 15 eV [1,12,37]. 

However, measurements from (In A) have given contradictory conclusions (see 
e.g. [38]): SAS [2] gives A (In A) = +1.0 between 2.5 - 6.3 PeV and constant 
above 6.3 PeV; HEGRA [39] gives a composition consistent with direct mea- 
surements at 100 TeV and gradually becomes heavier around ~ PeV; CASA- 
MIA [38] finds that (In A) increases in the range 100 TeV - 10 PeV; CASA- 
BLANCA [42] observes a composition becoming lighter between 1-3 PeV 
then heavier above 3 PeV. As pointed out in [8], these apparent discrepancies 
could as well be related to the different interaction models used (see Tab. 1 
and Fig. 7 of [42] where (In A) is displayed for four interaction models). At 
present QGSJET and VENUS seem to be favored [42], but it is clear that these 
Monte Carlo simulations are crucial to extract observables [8,43-45]. 
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Table 1 

Central values given by several experiments using several Monte Carlo simulations 
for hadronic interactions. 



Experiment [ref] 


Simulation 


-^knee 


7i 


72 


HEGRA [39] 


QGSJET 


4.0 PeV 


2.72 


3.22 


Tibet array [40] 


GENAS 


~ 2 PeV 


2.60 


3.00 


EAS-TOP [41] 


HDPM 


3 PeV 


2.76 


3.19 


CASA-BLANCA [42] 


QGSJET 


2-3 PeV 


2.72 


2.95 


KASCADE [8] 


QGSJET 


5.5 PeV 


2.77 


3.11 


KASCADE [8] 


VENUS 


4.5 PeV 


2.87 


3.25 



To sum up, from direct experiments and ground arrays, we can confidently 
present the following results: there is no break in spectra before a few hundreds 
of TeV [1] where usual cosmic ray acceleration is at work. Above a few PeV, 
the all-particle spectrum asymptotically reaches the slope ~ 3.0 (see Tab. 1) 
extending up to a second break [46]. Data show a gradual steepening of the 
spectrum rather than a single kink, but still the steepening happens within 
about one decade of energy. Focusing on the average logarithmic mass, it is 
found that (In A) is about constant near the knee and then gradually increases 
above the break (a change from a heavy to a light composition is then observed 
in the energy region 5 x 10 18 eV giving support to a different origin for these 
cosmic rays, i.e. extragalactic [46]). 

Actually, Stronger constraints can be obtained by combining the all-particle 
spectrum and the average mass: in [6], it is found that a simple model with a 
universal slope and a break at a given rigidity can match either the all-particle 
spectrum or the average composition, but not both. Among four models tested, 
[2] find that only one model (an adaptation of Biermann's one) is able to fit 
both observables at the same time. 

4-1.3 Tentative model for the knee 

At this stage, the best way to explain the data is to generate this break 
at a given energy or rigidity or energy per nucleon with a change in slope 
A7 ~ —0.4. Such a parameterization is sufficient to see what happens in our 
diffusion model. From what concerns the possibility of a different slope for 
p and other nuclei below the knee, we will set a p — a^>\ = 0.05, althought 
greater values are possible (aij denotes source spectra, whereas 7^ corresponds 
to propagated/measured spectra). Of course, the way we generate the break 
is to some extent unphysical. It should be kept in mind that the propagated 
all-particle spectrum and average logarithmic mass - which already present 
fine structures in this idealized case {i.e. sharp source spectral index break, 
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see results in fig7 ) would be even more complicated in a more realistic model. 
Anyway, the exact form of this break is certainly very model-dependent, even 
for the particular case of transition between different accelerating sources, 
because of the variety of source that can be invoked. 

4-2 Prescription for the propagation parameters 

General solutions and discussions about diffusion models can be found in [47]. 
The model we use is cylindrically symmetric with two zones - thin disc and 
diffusive halo - described at length in [14] and also used in subsequent pa- 
pers for primary [48] and secondary antiprotons [49], and for /^-radioactive 
nuclei [50]. The five parameters of this diffusion model are diffusive halo scale 
height L, rate of reacceleration related to the alfvenic speed V a , galactic wind 
V c , normalization Kq and slope 5 of the diffusion coefficient (see Sec. 1.1). A 
consistent range for the five parameters of our diffusion model has been ob- 
tained [14], but there is a strong degeneracy between the various parameters 
derived. Moreover, these parameters do not reproduce very well primary pro- 
ton fluxes. A solution for this problem is seeked (Donato et al, in preparation) 
adding a sixth free parameter to the study, namely the source spectral index 

(Xj. 

In this following, we focus on two specific parameters: the halo scale height L 
and the slope of the diffusion coefficient 5. The first one is related to geomet- 
rical effects since a thin halo (say ~ 3 kpc) would correspond to more "local" 
sources - in diffusion models, cosmic rays cannot come from regions whose 
distance is greater than that of the closer edge -. The second one is related to 
the source spectra a, since the measured spectrum slopes at Earth (7 ps 2.8) 
are linked via the approximate relation a ~ (2.8 — 5), 5 being the diffusion 
coefficient slope. To extend these calculations above the knee, we have two 
possibilities: (i) keep the same parameters throughout the energy range, and 
explain the knee by a change in slope spectra (see previous section), or (ii) 
explain the knee by a change of diffusive regime. We previously noticed that 
if the transition in diffusion regime is at a fixed rigidity, the situation is thus 
strictly equivalent to a break in spectra at E^ nee = Z eV (see Sec. 3). The 
transition proposed by Ptuskin and collaborators [51] is rather smooth, but 
as it requires 5 ~ 0.2 that is excluded by Maurin et al.'s analysis [14], we will 
not consider anymore this possibility. 

4-3 Radial distribution of sources q(r) 

Measurements of galactic 7 rays in the seventies have raised the question of the 
radial distribution of cosmic rays. This distribution is needed in order to eval- 
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uate the resulting gamma emissivity at different galactocentric locations. The 
first distribution used was that of Kodaira (1974) [52] following the radial dis- 
tribution of supernovae which is also close to that of pulsars. This is consistent 
with the present picture of cosmic rays where supernovae provide the energetic 
power and mechanism to accelerate nuclei. The description of the galactocen- 
tric distribution has been improved thanks to new observations of pulsars and 
supernovae. We take here the distribution of Case & Bhattacharya [53] which 
is an improvement of their earlier analysis [54]. Their last result provides a 
flatter distribution than the previous one, closer to the distribution adopted 
by Strong & Moskalenko [55]. A third form that we finally retain is a constant 
radial distribution: this will serve as a reference to estimate the pure radial 
distribution effect on the average logarithmic mass. These three distributions 
are reproduced in Tab. 4 and are presented in the left panel of Fig. 4 with the 
effects of metallicity gradient (see next section). 

4-4 Observed metallicity gradient 

The existence of radial metallicity gradients is now well established in spiral 
galaxies. Early studies showed a gradient for O/H from observations of ionized 
nebulae in galaxies like M33, M51, and M101, but later work observed this 
trend in our Galaxy for many other abundances (see [56] for a review and 
Sec. 2 of [57]). Several recent observations (see Tab. 1 of [58] for a compilation 
of results) lead to very similar conclusions for the metallicity gradient (see 
Table 2). 

Table 2 



Example of observed radial gradients (various elements) for several methods using 
samples distributed over the galactocentric distances vqc- 



Gradient 


Element 


Range 


Tool 


Ref. 


(dex kpc" 1 ) 




(r G c in kpc) 






-0.02 / -0.05 


C, O,... Gd 


6-11 


Cepheids 


[59] 


-0.07 (± 0.01) 


C, O, Mg, Si 


6-18 


Early B-type stars 


[60] 


-0.055 (±0.015) 


O, Ne, S, Ar 


4-10 


Planetary nebulae 


[61] 


-0.09 


Fe 


8-16 


Open clusters 


[62] 


-0.07 


N, O, S 


0-12 


HII regions 


[63] 



We choose the most recent study [59] which takes into account more than 25 
gradients corresponding to species from carbon to gadolinium (see in particular 
their Fig. 6-9 and Fig. 10). For all ions X (except He), the gradient is about 
d[X/K]/dr = -0.05 dex kpc" 1 . 
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4-5 Recapitulation of the model tested 



Let sum up the different models studied : 



(1) Source spectral index a - see Sec. 4.1: following our previous study [14], 
we set a = 2.8 — 5. 

Table 3 

Models I and II used in this paper to describe the slopes of cosmic ray sources, plus 
III and IV for the description of the knee. 



Model I 
Model II 
Model III 
Model IV 



All species with the same a 



H with a p , all others with a 



a below the knee and a% = a — 0.4 above the knee 
II below the knee, and a.2 = u — 0.4 for all species above the knee 



(2) Propagation parameters - see Sec. 4.2: three models with 5 (slope of the 
diffusion coefficient) equals to 0.46, 0.6 and 0.75 at a fixed diffusive halo 
size L — 10 kpc; one more with L = 3 (5 = 0.6). Exact numbers for 
other parameters are unimportant, however we emphasize that for each 
value of L and S, the parameters V a , V c and K are set such as to fit B/C 
data [14]. 

(3) Geometrical effects: 

• Radial distribution q(r) - see Sec. 4.3: three models. 

Table 4 

Models a,b and c used in this paper corresponding to three possible radial distribu- 
tion of sources. 



Model a 
Model b 

Model c 



q(r)= 1. 



W) ■ 
q(r) 



r \2.o 



.5 

r \ 0.5 



exp 



8^ 6XP 



-3.53 x 



-1. x 



(r-8.5) 

8^5 
(r - 8.5) 

8^5 



Case &: Bhattacharya [53] 
Strong & Moskalenko [55] 



• Metallicity gradient - see Sec. 4.4: as one can see in formula (11), 
there should be in the above formula an additional multiplicative factor 
(gx/Qp)© for ah species X. It is implicitly taken into account since all 

Table 5 

Two models - one with and one without gradient - used in this paper. 



Model 




No gradient 


Model 1 


Substitution q{r) - 


-> icr a05 ( r_8 -) x q(r), except for H and He 



fluxes are normalized to HEAO-3 data [64] at 10.6 GeV/nuc, except 
for p and He that are respectively normalized to ams proton data at 
79.6 GeV [65] and helium data at 49.2 GeV/nuc [66]. 
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Model Ia-0 will denote a model where all the sources have a fixed spectral 
index, the radial distribution is constant, and where there is no gradient. 
Model IIb-1 will denote a model where source spectrum index of H is different 
from all others, where q(r) is Case & Bhattacharya's one (see above), and 
where we choose a composition gradient of -0.05 dex kpc -1 for all species. We 
use this convention for all the figures proposed below. The reference model 
from where other will be compared is Model Ia-0 with 5 = 0.6, L = 10 kpc 
and no break. 



5 Results and concluding remarks 

All outputs come from the diffusion model of [14] with slight modifications of 
the inputs according to the above-prescriptions. 



5.1 Geometrical effects 

Curve Ia-0 (Fig. 4) corresponds to a model where all the species have the same 
spectral index, and where no geometrical effect is allowed. The first strong con- 
clusion is that the pure propagation effect affects dramatically the composition 
of cosmic rays. It is consistent with the conclusion of Sec. 2.2, namely that the 
evolution of the average mass is closely connected to propagation properties. 
Furthermore, at sufficiently high energy, as expected, we reach the asymptotic 
regime where propagation ceases to affect (In A). 

Other curves correspond to the evaluation of the pure geometrical effects. Both 
radial distribution (model a,b,c) and metallicity gradient (model 0,1) are sep- 
aratly presented. In fact, comparing curve Ia-0 and Ia-1, we could conclude 
prematuraly that metallicity effect plays a role in the evolution of (In A). How- 
ever, when the radial distribution of sources is correctly taken into account, 
metallicity only has a little additional effect (compare curves Ib-0 and Ib-1). 
Finally, the impact of metallicity is more or less pronounced depending of the 
distribution q(r) chosen. Nevertheless, as was correctly guessed in Sec. 2.3, (i) 
metallicity effect is of little importance and (ii) total geometrical effects cor- 
respond to an additional change of at most 5% compare to a non geometrical 
model. Furthermore, Fig. 4 corresponds to a halo scale height L = 10 kpc, for 
which the effects are maximized. For L = 3 kpc, these geometrical effects are 
completely negligible. 
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Radial distribution 



Average logarithmic mass 




Fig. 4. Left panel: representation of the three radial distribution of sources: solid lines 
are for constant distribution (models a), dotted lines are for Case & Batthacharia 
(models b) and dashed lines are for Strong & Moskalenko (models c); q(r) is evalu- 
ated (see Tab. 5) with (models -1, thin line) or without (models -0, thick line) metal- 
licity gradient. Right panel: average logarithmic mass (model I, 5 = 0.6, L = 10 
kpc) for the six cases presented on left panel. 



5.2 Effects of propagation and source spectra 



Fig. 5 (left panel) shows that the cosmic ray composition is very sensitive to 
the diffusion parameters, the strongest dependence being that of the diffusion 
power spectrum 5 (the fixed point around 500 GeV is just an artifact due 
to the normalization adopted). Note that the asymptotic diffusive regime is 
reached faster for larger values of 5. We checked that the influence of the 
parameter L is minor compared to <5's. We note that secondary contribution 
(nuclei produced by spallation of the main primary species) is important (right 
panel). In particular, when making the junction between direct measurements 
and ground arrays, these secondaries are almost never taken into account in 
the calculation of the average logarithmic mass, whereas they are implicitly 
counted in air showers data (their presence in the reconstructed quantity (In A) 
is more questionable [8]). This question is related to the ability of obtaining a 
confident normalization of (In A) with data from nuclear interaction models, 
e.g. [42]. 

The impact of the spectral difference between protons (slope cv p ) and other 
species (slope cvz>i = a) is illustrated in Fig. 6 (left panel) showing an increase 
of the mass composition with a p — a. The right panel summarizes the astro- 
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Fig. 5. Average logarithmic mass (model I, L = 10 kpc) for three values of the diffu- 
sion coefficient slope 5. Left panel displays these three values for two cases, (i) model 
Ia-0 (no geometrical effects, thin line) and (ii) model Ib-1 (Case & Bhattacharia's 
q(r) with V[Fe/H] = —0.05 dex kpc" 1 , thick line). Right panel shows for Model Ib-1 
the same three 6 values but introducing the distinction between (In A) evaluated (i) 
with primary species only (thick line), (ii) with all primary and secondary nuclei 
(thin line). 

physical effects studied here in the framework of a diffusion model: conclusions 
are similar to what was drawn in Sec. 2: pure propagation effects (Ia-0, dotted 
line) are mostly responsible for the increase of (In A) vs energy, geometrical ef- 
fects (Ib-1, dashed line) are less significant, and source effects (IIb-1, solid line) 
turn from an almost asymptotically constant (In A) into a constant enhance- 
ment of the same quantity. The three upper curves demonstrate importance 
of the diffusion power spectrum 5 and emphasize the role of the secondaries 
in the normalization of the average logarithmic mass. We also display the av- 
erage logarithmic mass as measured by experiments (stars). We renormalize 
to the observations at 100 GeV. We see that a difference a p — a > or/and 
large values of S are preferred. 



5. 3 Combination of all previous effects with a model for the knee 



In Fig. 7, we generate the knee either with a break at a fixed rigidity, i. e E^ nee = 
Z x 4 PeV, or at a fixed total energy per nucleus i. e E knee = A x 4 PeV. The 
possibility of a break at a single energy is not considered because it exhibits a 
very sharp break in all-particle spectrum not present in data. Fig. 7 displays 
the resulting curves for the two models IIIb-1 and IVb-1 that only differ in 
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Fig. 6. Left panel: average logarithmic mass (model IIb-1, 5 = 0.6, L = 10 kpc) 
for various values of cosmic ray proton spectral index a p (see Tab. 3). Right panel: 
summary of the various effects investigated in this paper for primary species, i.epure 
propagation effect (Ia-0) + geometrical effect (Ib-1) + spectrum effect (IIb-1). The 
three upper curves correspond to three 5 where secondary contributions are taken 
into account in Model IIb-1 (a p = a + 0.05). Stars are values of (In A) measured 
in direct experiments (taken from Fig. 29 of [1]) interpolated from their original 
normalization to our 100 GeV normalization value. 

their spectral indexes below the knee (a p — a = 0.05 for model IVb-1 and 
a p — a = 0. for IIIb-1). Compared to the results of the two components model 
(see Sec. 3), we remark both in (In A) and $ al1 some additionnal bifurcations 
generated by the helium component with a second even smoother transition 
provided by the CNO group. If these effects are not very relevant for the 
average mass composition, they hamper the interpretation of the transition 
from one regime to the other in the all-particle spectrum data. Secondaries 
smooth even more these transitions. They are also important for normalization 
of (In A). The two cases E kncc = Z x 4 PeV and E kncc = A x 4 PeV could 
be differentiated mostly through the all-particle spectrum. Finally, the bigger 
the difference between proton slope and other species before the knee, the 
smoother the bump in (In A) . 

Before concluding, we would like to make a brief comment on Swordy's model [67]: 
in the latter, an enhancement of light nuclei is predicted before the usual en- 
richment in heavier nuclei at the knee (some data can support this upturn, 
e.g. CASABLANCA's data, see Fig. 7). However, this requires a change in 
diffusion (i.e A5 < because $ oc R-(f= a + s )) at an energy smaller than the 
knee's, so that it should produce a bump visible in all-particle spectrum: the 
larger AS, the sharper the bump, and moreover a larger value of Aa is then 
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Fig. 7. Average logarithmic mass (left panel) and all-particle spectrum (right panel) 
for Model IIIb-1 and IVb-1 with a break either in rigidity (R = 4 PeV) or in total 
energy per nucleus (E/nuc = 4 PeV). Solid lines correspond to primaries only, 
whereas dotted lines correspond to primaries plus secondaries. In the right panel, 
fluxes for Model IVb-1 have been divided by two in order to avoid the overlap 
with Model IIIb-1. For illustrative purpose, some data from ground arrays have 
been displayed: triangles are CASABLANCA's data with Monte-Carlo HDPM [42], 
empty circles are from KASCADE collaboration with QGSJET simulation [8] (empty 
squares are from JACEE direct experiment; RUNJOB data plus some JACEE data 
are lower than 1.8 and do not appear on the graph, see [1]). 



necessary to reproduce the all-particle spectrum at the knee. Thus, if posi- 
tion and sharpness of this change in diffusive regime (quite constrained by 
all-particle spectra) is not theoritically excluded, it requires to be extremely 
fine tuned. Anyway, this model along with our theoretical predictions cannot 
be tested due to the large data scattering (see Fig. 7). 

The best clues up to now about the knee puzzle come from spectral analysis: for 
example the KASCADE collaboration [68] find that the all-particle spectrum 
exhibits a knee (A7 ~ 0.2 — 03) around 4 PeV, but that this knee is seen only 
in their light ion subsample for which A7 ~ 0.5. As regards the heavy ions, 
they find no changes in the region 1-10 PeV but the slope below the knee 
is smaller than that of light component. If this observation is confirmed, the 
average logarithmic mass is likely to evolve as depicted in our Model IVb-1. 
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5.4 Concluding remarks 



We have presented an analysis of the phenomena that affect the chemical cos- 
mic ray composition up to highest "Galactic" energies. Several astrophysical 
effects have been considered and geometrical effects have been found to play a 
minor role, while propagation effects (mostly selective destruction in flight of 
heavy nuclei) drive this evolution up to the knee where they cease to be effec- 
tive. A difference between the source spectrum of protons and other ions lead 
to a constant enhancement of (In A) up to the knee and ceases if the slopes 
above the knee are similar for all species. In the framework of a simple break 
in rigidity (or total energy per nucleon), a bump in the chemical composition 
occurs at the knee, but the larger the spectral difference between protons and 
other species, the smoother the bump. The secondary species induce an en- 
hancement in (In A) of about 15%. As a by-product, our study validates the 
approach recently used in [2] , i. e we demonstrated that above PeV energy, a 
propagation model including only sources and diffusion is relevant. Finally, 
the main problem of the diffusive problem is the normalization of the fluxes. 
It could be improved thanks to a better determination of the propagation 
parameters with more precise low energy data. 
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